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Abstract 


\ 


An  investigation  of  the  spin-orbit  resonance  effects  for  a  semi- 
synchronous  near  polar  orbit  was  undertaken  in  order  to  determine 
whether  consideration  of  all  resonance  terms  associated  with  a 
particular  commensurability  ratio  would  result  in  the  existence  of 
multiple  equilibrium  points  for  the  placement  of  satellites  utilizing 
this  type  of  orbit.  The  Hamiltonian  of  the  geopotential  in  Delaunay 
elements  using  first  and  second  order  zonal  harmonic  terms  (-#*-»—  J§,  aird 
J^^was  first  transformed  to  a  set  of  modified  variables.  The  effects 
of  the  resonant  disturbing  function  were  then  developed  and  the 
resultant  Hamiltonian,  valid  for  near  polar  inclinations  and  small 
eccentricities,  was  reduced  to  a  single  degree  of  freedom  through  a 
second  transformation.  Phase  portraits  of  the  system  Hamiltonian  were 
then  generated  and  an  analysis  of  the  resonance  structure  and 
librational  periods  performed. 

A  single  equilibrium  point  for  each  resonance  term  was  discovered. 
Because  of  the  physical  meaning  of  the  critical  argument,— the  stable 
equilibrium  points  of  both  the  first  and  second  resonance  terms  appear 
to  be  candidates  for  deployment  of  future  satellite  systems.  HEn-v' 
-«44it-iod>  the  semi-analytic  technique  of  phase  portraits  was  shown  to  be 
a  feasible  approach  to  the  investigation  of  resonance  effects. 
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I  Introduction 


The  advent  of  the  Space  Age  in  1957  brought  with  it  recognition  of 
the  inherent  advantages  of  observation  platforms  in  orbit  around  the 
earth.  Of  interest  to  meteorological,  land  resources,  and  navigation 
satellites  is  a  near  polar  orbit,  which  permits  coverage  of  all 
latitudes  of  the  earth's  surface  on  a  regular  basis.  It  has  also  been 
recognized  since  the  early  days  of  space  flight  that  satellites  are 
subject  to  significant  perturbations  due  to  the  non-spher ical  nature  of 
the  earth's  potential.  In  particular,  satellites  having  motions 
commensurate  with  the  rotation  rate  of  the  earth  experience  a  phenomenon 
known  as  spin-orbit  resonance  due  to  the  elliptical  shape  of  the  earth's 
equatorial  section.  Probably  the  best  known  and  most  studied  example  of 
this  phenomena  is  the  existence  of  two  stable  equilibrium  points  for 
satellites  in  geosynchronous  equatorial  orbits.  However,  resonance 
effects  have  also  been  observed  for  other  commensurate  orbits.  Of 
particular  interest  for  satellites  in  a  polar  inclination  is  a  semi- 
synchronous  or  12-hour  orbit.  As  the  missions  mentioned  above  generally 
require  the  utmost  in  orbital  stability  for  purposes  of  data  accuracy, 
any  stable  resonance  points  associated  with  this  general  type  of  orbit 
would  be  of  extreme  interest  for  positioning  future  satellites. 

Historical  Background 

A  review  of  the  literature  reveals  that  an  extensive  amount  of 
research  has  been  done  concerning  the  effect  of  the  earth's  shape  on 
satellite  orbits.  In  particular,  the  effect  of  the  earth's  oblateness, 
or  the  Ja  zonal  harmonic,  has  been  well  established  by  several 
researchers.  In  a  classic  edition,  the  November  1959  issue  of  The 
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Astronomical  Journal  contained  three  different  and  independent 
treatments  of  the  "main  problem"  of  satellite  theory. 

An  offshoot  of  these  early  studies  was  the  problem  of  motion  near 
the  critical  inclination  (63.4°)  which  can  be  considered  a  type  of 
resonant  motion.  This,  too,  received  much  attention  and  is  now  well 
unders  tood . 

Further  investigations  have  attempted  to  determine  the  effects  of 
the  earth's  longitude  dependent  harmonics  on  satellite  motion.  Because 
of  its  use  for  communications  satellites,  the  geosynchronous  orbit  has 
proven  to  be  of  the  greatest  interest. 

Blitzer  (Ref  7,8),  using  a  linearized  theory,  has  discussed 
equilibrium  solutions  for  a  geosynchronous  satellite  in  a  circular, 
equatorial  orbit  under  the  influence  of  the  principle  longitude 
dependent  term,  J2J.  Blitzer  (Ref  4,5)  later  treated  the  effect,  due  to 
higher  order  tesseral  harmonics,  on  a  geosynchronous  satellite  of  small 
eccentricity  and  inclination.  He  showed  that,  due  to  the  dominance  of 
the  J2l  term,  only  a  slight  displacement  of  the  equilibrium  positions 
occurred,  although  there  was  a  significant  change  in  the  librational 
periods. 

Musen  and  Bailie  (Ref  18),  in  1962,  studied  the  motion  of 
synchronous  satellites  incorporating  only  the  Jti  tesseral  harmonic  and 
the  Jt  and  J*  zonal  harmonics.  Their  analytic  technique  agrees  with 
Blitzer  for  synchronous  motion  in  the  equatorial  plane. 

Allan  (Ref  1),  in  1965,  discussed  the  motion  of  nearly  circular  but 
inclined  synchronous  orbits.  In  1967  (Ref  2)  he  studied  the  effect  of 
resonance  in  inclination  for  synchronous  satellites  in  nearly  circular 
orbits  when  the  positions  of  the  nodes  repeat  relative  to  the  rotating 
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primary.  Later  that  same  year,  he  also  studied  the  effect  of  resonance 
in  eccentricity  and  inclination  for  a  synchronous  satellite  in  a  nearly 
equatorial,  eccentric  orbit  which  occurs  when  the  longitude  of  the  line 
of  apsides  repeats  relative  to  the  primary  (Ref  3). 

In  a  book  published  in  1966,  Kaula  (Ref  16)  developed  expressions 
for  the  resonant  disturbing  function  of  the  geopotential  in  terms  of 
inclination  and  eccentricity  functions  and  derived  general  expressions 
for  the  variation  of  orbital  elements  due  to  arbitrary  zonal  or  tesseral 
harmonics . 

In  a  series  of  articles  over  the  years,  Garfinkel  (Ref  13)  has 
developed  a  formulation  known  as  The  Ideal  Resonance  Problem  which  has 
found  wide  application  in  resonance  theory  by  treating  specific  cases  as 
perturbations  of  the  same. 

Far  more  limited  is  the  literature  dealing  with  12-hour  resonant 
orbits.  A  1961  article  by  Cook  (Ref  10)  deals  with  resonance  effects 
for  commensurate  orbits  in  general  but  is  restricted  to  orbits  of  low 
eccentricity  and  inclination.  Blitzer  (Ref  6),  in  1963,  gives  a  cursory 
treatment  of  12  and  36-hour  orbits  in  addition  to  the  synchronous  case 
but  again  only  for  the  circular,  near  equatorial  case.  Allan,  in  the 
previously  mentioned  studies  of  1967,  also  considers  the  problem  of 
commensurate  orbits  other  than  the  24-hour  case  by  calculating  some 
resonance  strengths,  librational  periods,  and  optimum  inclinations  using 
inclination  functions. 

A  1972  article  by  Wagner  (Ref  21)  studies  the  longitude  drift 
regimes  of  eccentric  12-hour  orbits  utilizing  the  dominant  terms  of  the 
geopotential.  Inclinations  near  the  critical  inclination  are 
specifically  examined  because  of  their  possible  long  term  stability. 
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Ten  years  later,  Sochilina  (Ref  20)  investigated  specific  cases  of 
12-hour  orbits  at  the  critical  inclination.  In  particular,  he  looks  at 
orbits  of  the  "Molniya"  or  "Navstar"  type. 

From  the  preceding,  it  is  apparent  that,  while  considerable 
progress  has  been  made  in  resonance  theory,  it  appears  at  this  time  that 
little  hope  exists  for  a  general  analytic  solution  except  for  certain 
specific  cases.  These  include  mostly  satellites  whose  mean  motion  is 
strictly  commensurate  with  the  earth's  rotation  rate  a*nd  whose  orbits 
are  at  the  critical  inclination  or  have  zero  eccentricity.  It  should  be 
noted  that  while  Dallas  and  Diehl  (Ref  11)  claimed  to  have  found  such  a 
general  solution  in  1977,  they  were  later  shown  by  Jupp  (Ref  15)  to  have 
made  a  serious  error.  Garfinkel  (Ref  13),  in  his  1979  summary  of  the 
Ideal  Resonance  Problem,  lists  the  synchronous  satellite  with  non-zero 
eccentricity  and  the  general  p:q  resonance  between  the  period  of 
revolution  of  the  satellite  and  the  rotation  of  the  earth  as  two  of  the 
outstanding  unsolved  problems  of  resonance  theory. 

As  a  result,  although  the  effort  continues  to  find  formal,  global 
solutions  to  resonance  problems,  most  recent  efforts  address  themselves 
to  specific  satellite  orbits  such  as  the  studies  by  Wagner  or  Sochilina 
or  attack  the  problem  via  a  semi-analytic  or  numerical  technique  such  as 
used  by  Nacozy  and  Diehl  (Ref  19)  in  a  recent  article. 

Of  this  latter  category  is  the  method  of  computer  plotting  phase 
portraits  of  the  specific  system  Hamiltonian.  Although  this  has  been 
done  previously  in  order  to  provide  some  sort  of  "physical  feel,"  it 
would  appear  that  little  has  been  done  to  utilize  the  process  as  the 
primary  investigative  tool.  While  this  technique  provides  information 
only  in  the  regime  of  the  particular  inclination  studied,  it  is  valid 
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for  small  eccentricities,  permitting  investigation  of  all  resonance 
terms  associated  with  a  particular  commensurabi 1 ity  ratio  and  is  easily 
modified  to  incorporate  additional  harmonic  terms  or  study  different 
commensurabilities.  But  most  importantly,  it  circumvents  most  of  the 
mathematically  sophisticated  and  algebraically  laborious  techniques 
required  in  the  search  for  more  general  analytic  solutions. 

Objective 

In  an  effort  to  better  understand  the  regime  of  the  semi- 
synchronous  polar  orbit,  this  study  investigates  the  resonance  effects 
due  to  the  combination  of  certain  zonal  and  tesseral  harmonics  of  the 
geopotential  via  computer  generation  and  analysis  of  phase  portraits  of 
the  sys tem  Hamiltonian.  It  was  suspected  that  consideration  of  all 
resonance  terms  associated  with  a  particular  harmonic  term  and 
commensurabil ity  ratio  as  a  result  of  non-zero  eccentricity  would 
result  in  multiple  equilibrium  points.  And,  although  not  truly  semi- 
synchronous,  these  points  could  sufficiently  reduce  station-keeping 
requirements  so  as  to  merit  careful  consideration  for  the  deployment  of 
future  satellite  systems. 

Scope 

Although  the  technique  is  easily  modified  to  incorporate  additional 
tesseral  harmonics  or  commensurabil ity  ratios,  the  scope  of  this 
investigation  is  limited  to  the  resonance  effects  arising  from  the 
interaction  of  the  first  and  second  order  zonal  effects  of  the 
geopotential,  i.e.,  the  J2,  J|,  and  J\  terms,  with  the  principal 
longitude  dependent  term,  J2i,  on  a  semi-synchronous,  near  polar  orbit. 
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General  Approach 


The  method  of  investigation  revolves  around  utilizing  a  two-body 
Hamiltonian  of  the  system  that  already  contains  the  zonal  harmonic 
terms.  The  Hamiltonian,  which  is  expressed  in  terms  of  Delaunay 
variables,  is  then  converted  by  a  canonical  transformation  to  a  set  of 
variables  that  will  prove  more  convenient  in  the  development  of  the 
disturbing  function.  The  appropriate  resonance  terms  of  the  J22 
tesseral  harmonic  are  then  selected  and  transformed,  assuming  small 
eccentricity,  for  inclusion  in  the  system  Hamiltonian.  Next,  a  second 
canonical  transformation  is  used  to  eliminate  variables  and  reduce  to  a 
single  degree  of  freedom.  The  computer  is  then  used  to  plot  phase 
portraits  of  the  Hamiltonian  of  each  resonance  term.  These  show  the 
number  and  structure  of  the  stable  equilibrium  points.  In  addition,  the 
program  determines  the  width  of  the  stable  region,  the  libration  period 
about  the  stable  point,  and  its  location  relative  to  a  nominal  semi- 
synchronous  orbit. 

Sequence  of  Presentation 

Chapter  II  presents  the  theoretical  development  of  the  system 
Hamiltonian  and  equations  of  motion  for  each  resonance  term.  Chapter 
III  discusses  the  computer  analysis,  Chapter  IV  presents  the  results, 
and  Chapter  V  provides  an  interpretation  and  a  brief  discussion  of  the 
results  as  well  as  recommendations  for  further  study. 
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II  Theoretical  Development 


Phase  Portraits 

Prior  to  actually  developing  the  equations  necessary  for  this 
investigation,  some  'justification  is  in  order  by  way  of  a  brief 
discussion  of  the  theory  of  phase  portraits. 

In  Hamiltonian  mechanics,  the  motion  of  a  dynamical  system  with  n 
degrees  of  freedom  can  be  represented  geometrically  by  the  motion  of  a 
representative  point  P  in  a  2n-space  defined  by  the  n  generalized 
coordinates  and  the  n  conjugate  momenta  p^.  This  2n  space,  called 
the  phase  space,  can  be  considered  a  2n  dimensional  Euclidean  space. 

The  physical  difference  between  q^  and  p^  poses  no  problem  since  they 
possess  equal  status  in  Hamiltonian  mechanics. 

In  the  phase  space,  the  point  P  traces  a  path  or  trajectory  that 
corresponds  to  a  particular  solution  of  the  dynamical  problem.  For 
different  initial  conditions,  the  motion  will  be  described  by  a 
different  trajectory.  The  real  advantage  of  the  geometric 
representation  becomes  evident  when  we  consider  not  just  one  trajectory 
but  the  totality  of  trajectories,  known  as  the  phase  portrait, 
representing  all  possible  solutions. 

Of  course,  the  representation  of  the  motion  of  a  system  with  n 
degrees  of  freedom  in  a  2n  dimensional  phase  space  is  itself  difficult 
and  hard  to  visualize.  If  the  system  could  be  reduced  to  a  single 
degree  of  freedom,  the  motion  could  be  portrayed  on  a  two  dimensional 
phase  plane  which  is  easier  to  comprehend.  Hence,  if  constants  of  the 
motion  could  be  obtained  by  elimination  of  variables,  a  complex 
dynamical  system  could  be  easily  portrayed  on  a  two  dimensional  phase 
portrait.  This  can  be  done  relatively  easily  by  performing  a  canonical 
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transformation  of  variables. 

The  goal  of  much  of  the  rest  of  this  chapter  is  to  develop  the 
Hamiltonian  for  a  semi-synchronous  near  polar  orbit  and  by  a  canonical 
transformation,  reduce  it  to  a  single  degree  of  freedom  so  that  the 
resonance  effects  can  be  portrayed  on  a  two  dimensional  phase  portrait. 


Zonal  Hamiltonian 

In  order  to  study  the  resonance  effects  for  any  type  of  orbit,  the 
Hamiltonian  for  a  two-body  system  that  already  contains  the  effects  of 
the  included  zonal  harmonics  must  be  obtained.  Ignoring  the  long  period 
terms,  Hori  (Ref  14:292),  after  the  method  of  Brouwer  (Ref  9),  gives  the 
Hamiltonian  through  second  order  as 

F  =  F,  +  F,  +  F*  (1) 

with 


(2) 


F,  = 


yHj’Re  1  3  H 

T75T<-  7  ♦  f<§)2) 


2L’GS 


(3) 


1** JiR*  15  27  .27  135 

4LSG* ' 32  32  J*  “  U6  16  J|UG; 


,15  .  315 


v32  32  Jf;V  8  G'*  ”'G'  '  "  'G 

.  (ij.rii  +  «  J_1  _  (15  +  225  _  ( 105  _  525  £\(H)M> 

V  32  32  J|  U6  16  Jf'V  '  32  32  J|JV  J> 


1  kl  -  6(|)J  +  9(|)') 


(4) 


where  y  is  the  gravitational  parameter  of  the  earth,  its  mean 
equatorial  radius,  and  L,  G,  and  H  are  the  momenta  terms  of  the  Delaunay 
action-angle  variables.  These  are  related. to  the  standard  Keplerian 
elements  by 
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t  =  M  =  mean  anomaly, 


G  =  L/l  -  e 2 ,  g  =  to  =  argument  of  perigee, 

H  =  G  cosl,  h  =  0  =  right  ascension  of  ascending  node, 

where  a  and  e  are  the  osculating  semi-major  axis  and  eccentricity, 
respectively,  and  I  is  the  instantaneous  inclination.  The  subscripts  of 
F  represent  the  order  of  the  associated  zonal  harmonic,  J2  being  first 
order  and  J \  and  J„  of  the  second  order. 


First  Transformation 

It  will  prove  convenient  in  the  development  of  the  disturbing 
function  to  have  an  eccentricity  related  variable  that  is  itself  small 
for  small  eccentricities.  This  can  be  done  by  performing  a  canonical 
transformation  to  a  new  set  of  variables.  Defining  the  new  momenta 
variables  as 

X  =  L  =  ATa 

Y  =  L  -  G  =  /^7[1  -  /I  -  el] 

Z  =  H  =  (X  -  Y)  cosl 
the  generating  function  becomes 

S2  =  XI  +  (X  -  Y)g  +  Zh  (5) 

and 

x=t+g,  y=-g.  z=h 
become  the  new  angle  variables.  Converting  to  canonical  units 

(R  =  y  =  1),  the  zonal  Hamiltonian  now  becomes 

e 

F'  -  F' (X,Y,X,x,y,z, t)  =  F(L,G,H,l,g,h,t)  =  F,  +  F,  +  F*  (6) 

where 


F. 


(7) 


9 


(8) 


J  2  1  3  Z 

F|  ~  2X 5 (X-Y ) 3  ^ ~  ^  +  ^ 


2  X-Y 


F  =  --  —  /II  +  I?  _  ,27  135  ^)(_Z_)a 

*  4X 5  ( X-Y ) 5  ^  3  2  32  J’  V 16  16  J*nX-Y; 


+  (11  +  315 

V32  32  J2 


32  Jf  X-Y 


)(^7>'  +  - 


2 

3  X 


8  X-Y(1  ”  6(X-Y)J  +  9(X-Y)I,) 


+  «  -  (22  +  225 

X-y'  1 32  32  J|  16 


_  \  (  2  . 2 

16  jJnX-YJ 


,105  _  525  J\f_Z_x,n 
32  32  jT  X-Y 


(9) 


Disturbing  Function 

At  this  point,  the  resonance  disturbing  function  can  be  developed 
and  incorporated  into  the  Hamiltonian. 

The  disturbance  at  a  point  due  to  the  principle  tesseral  harmonic 
term  of  the  geopotential,  J2J)  is  normally  given  as 

w  Re 

U22  =  ^2 2^1  (s in$)cos2( X  -  X22)  (10) 

where  r,  8,  and  X  are  the  geocentric  radius,  latitude  and  longitude  of 

the  point,  X22  is  the  longitude  of  the  major  axis  of  the  earth's 

equator,  and  Pf(sinB)  is  the  associated  Legendre  function.  Inclusion  in 

the  Hamiltonian  usually  entails  expanding  Eq  (10)  in  powers  of  the 

eccentricity.  At  this  point,  however,  a  slight  short-cut  is  in  order. 

From  Allan  (Ref  1:62),  the  general  form  of  the  disturbing  function  in 

canonical  units  (R  =  el)  can  also  be  expressed  as 
e 

J  l 

U,  *  ■  -1”,-  l  F,  (I)G  (e)cos[  (t-2p)w  +  (l-2p+q)M  +  m(fl-n,t-X  )]  (11) 
Im  t+1  L  Imp  tpq  Im 

a  p,q  r 

where  a,  e,  I,  fl,  to,  and  M  are  the  Keplerian  elements,  n#  is  the  mean 
rotation  rate  of  the  earth  such  that  n,t  is  the  Greenwich  sidereal  time, 
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and  F,  (I)  and  G.  (e)  are  inclination  and  eccentricity  functions  for 
Imp  Ipq 

which  formulas  or  existing  tables  are  available  (Ref  16). 

To  pick  out  the  resonant  terms,  note  that  the  rate  of  change  of  the 
argument  of  Eq  (11)  is 

(fc-2p)u  +  (t-2p+q)(n+x  )  +  m(fi-n0) 

* 

where  X  is  the  modified  mean  anomaly  at  epoch  defined  by 

M  =  Jndt  +  x 

•  «  •  ,  »  • 
and  n  is  the  orbital  mean  motion.  Since  tu,  x  >  and  SI  are  small  for  a 

semi-synchronous  polar  orbit,  this  is  approximately  equal  to 

(l-2p+q)n  -  mn0 

For  resonance  to  occur,  this  must  be  equal  to  zero  and  since  n  s  2n0  for 
a  semi-synchronous  orbit,  the  resonance  condition  is  given  when 
m  =  2(l-2p+q).  Discarding  the  non-resonant  terms  leaves  the  following 
general  expression  for  the  resonant  terms: 

J  t 

U,  =  l  F.  (I)G  (e)cos  [  (l-2p)u>  +  (l-2p+q)M  +  mE]  (12) 

ixa  i+i  imp  ipq 
a  p 

where  E  =  fl-njt-A^  and  q  =  m/2-i+2p.  Therefore,  for  the  J2  2  term 
(l  =  m  =  2)  p  =  0,1,2  and  q  =  -1,1,3  which  gives  F2I#,  F221,  and  Fl22  as 
the  inclination  functions  and  G2,_,,  G2ll,  and  G22,  as  the  associated 
eccentricity  functions. 

From  Kaula  (Ref  16:34-38)  the  inclination  functions  are 

15  f(l  +  cosl)1  ,  F211  =  |sin2I  ,  F2  2  2  «  |(1  -  cosl)1 


and  the  eccentricity  functions  to  0(e7)  are 


G 


*»-i 


e  e1  5e}  I43e 7 . 
2  16  "  384  “  1'8432; 
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t 


„  ,3e  .  27e5  261es  143097e7 

G*“  '2  16  +  128  6144 


G 


2  2  3 


,e*  lie5  31 3e 

''48  768  30720 


7 


gives  for  the  resonant  disturbing  function  for  J 

3d «  )  3  c  S 

*  T.2/  e  e  5e5 

4P“(1  +  cosI)  (“  2  +  16 - 


2  2 


Se5  143e7 

384  -  T8432)C0S(M  +  2E  + 

3J*2  •  ,w3e  .  27e 5  .  261es 

iP“sln  I(r  +  ir  + 


3Jj; 


27e 5  .  261es  .  143097e\ 
v2_  ^  'T6_  “T28-  +  614~)c°s(M  +  2E) 


(13) 


It  is  worth  noting  at  this  point,  that  resonance  effects  on  a  semi- 
synchronous  orbit  only  exist  for  non-circular  cases  since  there  are  no 
zero  order  terms  in  the  eccentricity  functions. 

Before  adding  to  the  zonal  Hamiltonian,  the  disturbing  function 
must  now  be  converted  to  the  modified  variables.  Since 

Z  *=  (X  -  Y )cosI  (14) 

then 


cosl  =  ^  (15) 

and  letting 

B  -  ~  «  cosl  (16) 

gives 

(1  +  cosl)1  -  (1  +  B)1  (17) 

sin*I  =  1  -  B*  (18) 

(1  -  cosl)1  -  (1  -  B)*  (19) 

and  if 


I 
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(20) 


then 

e  =  /[2Y/X  ~  (Y/X)2]  (21) 

Now,  because  of  the  earlier  change  of  variables,  if  e  is  constrained  to 
be  small  means  that  Y  is  also  small  and,  as  a  result,  the  higher  orders 
of  Y/X  can  be  discarded  in  Eq  (21)  as  well  as  the  higher  orders  of  e  in 
the  resonant  disturbing  function. 

Making  the  appropriate  substitution  of  variables  into  Eq  (13)  gives 


2  1  _ 

U,,  -  ££T~(1  +  B)  2  (-  -|/2Y/X)cos(x  -  y  +  2z  -  2n,t  -  2A22) 


3J , 


+  2xT-(l  -  B2)(|/2Y/X)cos(x  +  y  +  2z  -  2n,t  -  2A22) 

3d  2  2  Y  i 

+  4Xr-(1  "  B)1(24x'/2Y/X)c0s(x  +  3y  +  2z  -  2n»t  ~  2*2i> 


(22) 


Simplifying  results  in 
3J22/y 

Ujj  = - - ( 1  +  B)2cos(x  -  y  +  2z  -  2n#t  -  2X22) 

4/2X1 * 

\ 

9J22/Y 

+  — ;  -  ■■■— ( 1  -  B2)cos(x  +  y  +  2z  -  2n,t  -  2X22) 

2/2Xrr 

3J2  2/y* 

+  - - (1  -  B)2cos(x  +  3y  +  2z  -  2n,t  -  2X21)  (23) 

48/2X 1  * 


Second  Transformation 


To  investigate  the  individual  resonance  terms,  each  one  is 
incorporated  separately  with  the  zonal  Hamiltonian.  This  is  possible 
because  for  the  polar  orbit  the  effects  of  the  individual  terms  are 
distinctly  different.  In  the  immediate  vicinity  of  each  resonance,  the 
effects  of  that  term  dominate  the  others  and  may  be  treated 
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independently.  In  order  to  do  this,  the  explicit  time  dependence  must 
be  removed.  This  can  be  accomplished  by  means  of  a  second  canonical 
transformation.  At  the  same  time  this  process  can  be  used  to  reduce  the 
Hamiltonian  to  a  single  degree  of  freedom  by  eliminating  variables. 

This  is  done  by  grouping  most  of  the  angle  variables  in  the  cosine 
function  together  into  a  single  critical  argument. 

In  order  to  make  the  transformation  perfectly  general  for  all  three 
resonance  terms  as  well  as  those  of  other  conrmensur abi lity  ratios  let 

s  =  ix  +  jy  +  kz  -  mn0t 
where  i,  j,  k,  and  m  are  integers, 

q  =  x 

and 


r  =  z 

The  generating  function,  S2,  then  becomes 

S2=Qx+S(ix+jy+kz-mnPt)  +Rz 


yielding 


X  =  Q  +  iS ,  Y  =  jS ,  Z  =  kS  +  R 
as  the  relationship  between  the  old  and  new  momenta.  The  new 
Hamiltonian  then  becomes 


with 


* 

F 


*  as2 

F  (Q,R,S,_,_,s,_)  =  F'(X,Y,Z,x,y,z,t)  - 
F,  +  Fi  +  F2  +  U22  +  mn,S 


(24) 


(25) 


*  2(Q+iS)‘ 


(26) 


F, 


Ji 

2 (Q+ iS  )  * 


A’ (-  |  ♦  |b*) 


(27) 
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15  27  J|'  27  135  J* 

J  4(Q+iS) 1 0  ' 32  32  Jf  l16  16  Jf' 


J, 


J  » 


*  <31  *  T§  TT>b*  *  tA(1  -  68‘  *  9B"> 


-  ...»  ,  il  il  _  ,11  .  225 

A  1 32  n  t  2  '■it 


16  +  "t f  J7^2 


32 

/105  525  J\rM% 
(  32  32  JI)B  1 } 


(28) 


where 


A  = 


Q+iS 


Q+(i-j)S 


B  = 


kS+R 


Q+(i-j)S 


and  Uj j  is  the  corresponding  resonance  term.  In  the  current  case,  only 
the  index  j  changes.  As  a  result,  the  only  changes  for  each  resonant 
Hamiltonian  are  the  actual  resonance  term  and  the  coefficients  A  and  B. 

Now,  since  q  and  r  do  not  appear  in  the  new  Hamiltonian,  Q  and  R 
are  constants  of  the  motion.  In  addition,  since  time  does  not  appear  in 
the  Hamiltonian,  the  total  energy  of  the  system  is  a  constant  also  and 
is  equal  to  the  Hamiltonian.  Hence,  the  system  has  been  reduced  to  a 
single  degree  of  freedom  where  phase  portraits  of  the  coordinate  s  and 
the  conjugate  momentum  S  for  various  values  of  the  Hamiltonian  can  be 
plotted  by  the  computer.  To  do  this,  however,  requires  determining  the 
location  of  the  resonance  equilibrium  points.  Since  the  equilibrium 
points  are  relative  extremas  of  the  system  Hamiltonian,  i.e.,  an  energy 
well  or  peak,  their  location  is  determined  by  where  the  derivatives  of 
the  Hamiltonian  are  equal  to  zero.  As  the  Hamiltonian  is  a  function  of 
both  S  and  s,  the  derivative  with  respect  to  both  of  them  will  be 
required.  These  derivatives  will  also  prove  necessary  in  the  actual 
computer  plotting  since  they  define  the  slope  of  the  Hamiltonian  at  a 
given  point. 
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Derivatives 


Note  that  since  s  only  appears  as  an  argument  of  the  cosine, 

* 


3F 

3s 


=  "  U, 


(29) 


with  the  cosine  replaced  by  the  sine.  The  other  partial  is  not  quite  as 

.  * 
straightforward.  In  order  to  save  some  effort,  3F  /3S  for  the  zonal 

part  can  also  be  done  in  a  general  form.  Liberal  application  of  the 

chain  rule  to  Eqs  (26)-(28)  and  a  lot  of  "attempted  simplification" 

results  in  the  following  expressions  for  the  partials  of  the  zonal  terms 


3F  o 


3S  (Q+iS)' 


(30) 


ill  =  -  3J*  .»fr(2i-j)Q*(2i«-2ij)Su  l  3  2) 

as  2(Q+iS)’  A  Q+( i- j )S ](  2  2®  ' 

_  AitrkQ~(i~j)Rn 

AB  Q+(i-j )S  ]> 


(31) 


5J* 


12.  _  _  A1  [  f  ( 2 i~ j  )Q+ (2i *-2i j )S.  .15  +  27  ^  27  +  U5  ^ 

3S  4(Q+iS)“  11  Q+(i-j)S  M32  32  J*  U6  16  J*'B 

♦  (H  *  i>8'  *  lA<1  -  6B'  *’*'>  -  -  if  3T 


.15  .  225 


2 

J* 


.105  525 


J «. 


*  Ia<-  128  +  36b»)  ♦  a«[(-4  +  ™  jf)B  +  (ir  -  ¥  ir)Bli> 


♦  [ 


_  1<2(,  _  6B2  +  9B*)  .  2Afil  +  45  _  ,15  225  f 

Q+(i-j)SJ<8U  6B  9B  '  2A  32  32  (16  16  J|)B 


J2_ 


-  (121  -  525 

(  32  32  jf;B  1>}1 


(32) 


The  specific  Hamiltonian  and  derivatives  for  investigation  of  each 
individual  resonance  term  can  now  be  determined. 
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First  Resonance  Term 


Taking  the  argument  of  the  first  term  of  Eq  (21)  gives 
transformation  indices  of  i  =  1,  j  =  -1,  k  =  2,  and  m  =  2.  Transforming 
the  resonant  term  and  applying  the  indices  to  Eq  (27)  results  in  the 
following  resonant  Hamiltonian 


* 

F  = 


1  *  J*.  A»,  1_  +3_rM 

2(Q+SP  2(Q+S)‘  A  2  2  B 

J'  *  17  J-  ,27  135  J\Rl 

4 (Q+S ) 1 4  A  {32  32  Jj  ( 16  16  J| >B 

*  'M  *  tI  jt)b'  *  lA<1  - 6B’  * 9B') 

.ir15  .  45  ,15  .  225  * 

_  A  [32  +  32  Jf  ‘  (T6  +  ~J6  jf)B 

/ 105  525  J* 

(_32  “  “32  JI)B  3 } 

- (1  +  b) 1cos(s  -  2X,j)  +  2n0S 

4/2  (Q+S )*» 


(33) 


=  Q*S 
Q+2S 


and 


_  2S+R 
Q+2S 


Likewise,  taking  the  partial  of  the  resonant  term  and  applying  the 
indices  to  Eq  (29)  gives 


li¬ 

as 


w  -  i  *  !»■)  -  ab[2^§]> 

K^yrr  *  §  £  -  <g  *  HI  > 

*  <M  *  m  ji)B-  *  I**1  - *•■  * ,B,)  -  *‘iM  *  35  jt 


(11  +  125  J'  t 
16  16  J|)B 


fK(5  _  525  1\bM> 
(  32  32  J*,B  J> 
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A/r2Q-2R..  ,27  135  J*  ,15  315  J*. 

~  5{I-fes]<'  ("8  +  —  II)B  +  (“8  +  —  If)BS 

♦  |a(-  12B  +  36BS)  +  A2[(^|  +  yr)B  +  (■—  _  ji)B‘]> 

"  [QT2SKI(1  "  6B*  +  9B%)  '  2A[if  +  M  Jf  ‘  (lf  + 

T. 


16  jf 


—  -  11 - [ ( 12S-Q)( 1  +  B)2  -  (2S~9o HSA( 1  +  B)]cos(s  -  2AZJ) 

8/-2S(Q+S)‘ 5  Q  S 


+  2n0 


and 

*  3Jil/^S 

~  =  (1  +  B)2sin(s  -  2AJ2)  (35) 

8  4^2 (Q+S)1 * 

Note  that  the  appearance  of  /-if  is  taken  care  of  by  the  fact  that  the 
specific  form  of  the  transformation  results  in  S  itself  being  negative. 
This  has  some  implications  for  the  computer  implementation. 


Second  Resonance  Term 

Taking  the  second  term  of  Eq  (21)  gives  indices  of  i  =  1,  j  =  1, 
k  =  2,  and  m  =  2.  The  only  part  of  the  Hamiltonian  that  changes  is  the 
resonant  term 

9JaJ/s 

U22  *  ■■  —  -(1  -  B2  )cos(s  -  2Xj  2 )  (36) 

2/2(Q+S) 1 4 

and  the  coefficients 


Q+S 

Q 


and 


2S+R 


The  partials  of  F  are  now 
* 


8F 

as 


1 


(Q+S ) * 


3J, 

2<Q+S)7 


A’{(-  j  +  |b2)  -  2AB> 
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1  .s.,15  .  27  ^  ,27  135  1\_2 

4 (Q+S ) 1 1  j2  32  J*  ~  U6  +  16  J| ' 

_  A.  .  (27  135  ^  (15  315.  J'  , 

5<2<  C-g  +  -g-  jr)B  +  (~g  +  ~g  -  jJ)B 

♦  |A(-  12B  ♦  36B1 )  +  AM  (if  ♦  2|5  ^  +  (105  _  5|5  ^)fil]> 
+  <|(1  _  6B,  «.  9B-)  _  2A(1|  +  ^|  ^  .  (||  +  125  ^)B, 


4/2S(Q+S) 


—  [  ( 12S-Q  K 1  -  B1)  +  8SAB]cos(s  -  2X22)  +  2n, 

1  5 


9J2  2/if 

- 1 

2/2CQ+S)15 


*(1  -  Bl)sin(s  -  2Xjj) 


Third  Resonance  Term 

The  third  term  of  Eq  (21)  gives  indices  of  i  =  1 ,  j  =  3 ,  k  =  2 ,  and 

* 

m  =  2.  Again,  the  only  change  in  F  is  the  resonant  term 


SJjj/ijs’ 

UM  =  - —  --(1  -  B)  *cos(s  -  2X2  2 ) 

48/2(Q+S ) 1 5 

and  the  coefficients 


A  -  -2+i 

Q-2S 


and  B  = 


The  partials  of  F  become 


-  W  -  5  *  !»■>  - 
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■  4 (q+s ) 1 1  aS[{_q4I><M  +  31  Jl  ■  (T6  +  jI)b2 

+  (M  +  ^  TF)b*  +  lA<1  -  682  +  5B%)  -  A’[M  +  31  71 


.15  .  225 


J. 


J* 


+  |A(-  12B  +  36B1 )  +  A*  t  (— ^  +  ^  jT)B  +  C1!1,-  ~  j|-)BJ]> 


.  r  3Q  3f  ,  g_.,  _  ,15  45  J*  _  ,15  225  J'  2 

Q-2S  <8  1  68  98  2A  32  32  71  16  16  7f)B 


,105  525 

(~32  "  “32  7T)B  ]>>1 


3J12/3S 
32/2(Q+S)i T 
+  2n0 


[  ( 12S-3Q)(  1  -  B) 2  +  (^||)4SA(1  -  B)]cos(s  -  2X„) 


(40) 


and 


3F 

3s 


9j2J/3s= 

48/2(Q+S) 1 


;(1  -  B)*sin(s  -  2X2i) 


(41) 


Equilibrium  Conditions 

Location  of  the  stable  and  unstable  equilibrium  points  is 

...  * 

determined  by  where  the  two  partial  derivatives  of  F  are  both  equal  to 
zero.-  By  inspection  of  Eqs  (35),  (38),  and  (41),  one  can  see  that  all 
three  resonances  will  have  equilibrium  points  where 

s  =  2Xj  j 

and 


8  *  w  +  2Xtl 

Determination  of  the  values  of  S  and  the  exact  structure  of  the 
resonances,  however,  will  require  computer  analysis  as  explained  in 
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Chapter  III.  Before  that  can  be  done,  though,  one  last  development 
needs  to  be  performed. 


Libration  Period 

Since  an  object  not  placed  right  at  a  stable  equilibrium  tends  to 
librate  or  oscillate  in  semi-major  axis  and  longitude  about  the  point, 
it  would  be  useful  to  know  the  libration  period  for  each  resonance 
region.  Since  only  the  region  around  the  stable  point  is  of  major 
interest,  only  the  period  for  small  displacements  near  the  stable  point 
need  be  calculated.  This  is  normally  done  analytically  by  use  of 
elliptic  integrals.  Again,  a  short  cut  is  available. 

In  the  neighborhood  of  an  equilibrium  point  of  a  function  such  as  a 
stable  point  which  is  a  minima  of  the  Hamiltonian,  the  equations  of 
motion  may  be  linearized  about  the  equilibrium  point  (Ref  17:180)  such 
that 

• 

6S  =  If 6s  (42) 

and 

6s  =  t§5S  (43) 


for  some  small  displacement  6s  or  6S  about  the  stable  point.  But,  from 
Hamilton's  equations 


§ 


(44) 


and 

* 

.  8F 

am  —  ■ 

as 

Substituting  Gqs  (44)  and  (45)  into  Eqs  (42)  and  (43)  yields 


(45) 
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Hence,  by  numerically  evaluating  the  second  partials  for  a  small 
displacement  near  the  stable  point  where  the  first  partial  is  known  to 
be  zero,  the  librational  period  for  small  oscillations  about  the  stable 
point  can  be  determined. 
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Ill  Computer  Analys is 

With  the  necessary  equations  now  developed,  the  next  step  was  to 
implement  the  computer  analysis  of  the  effects  of  the  three  different 
resonance  terms.  Because  of  the  extremely  complex  nature  of  the 
equations,  this  analysis  was  performed  on  a  large  CDC  mainframe 
computer. 

Sof tware 

The  appendix  contains  a  commented  listing  of  the  routines  (set  up 
for  the  first  term)  used  in  the  analysis.  They  were  implemented  in 
FORTRAN  V  because  of  the  availability  of  a  contouring  routine  in  that 
language  and  the  author's  familiarity  with  FORTRAN.  The  four  main 
routines  used  or  developed,  are  briefly  described  in  the  following 
sections. 

Program  RESPLT.  RESPLT  is  the  driver  program  for  the  software.  It 
is  responsible  for  run  initialization,  program  control,  and  phase 
portrait  plotting.  The  program  itself  can  actually  be  broken  down  into 
four  main  functions.  The  first  is  to  initialize  all  the  program 
constants  including  those  for  the  geopotential  model.  Next,  the  plot  is 
initialized  by  drawing  and  labeling  of  the  plot  axes.  Third,  the  stable 
and  unstable  equilibrium  points  as  well  as  resonance  width  and  libration 
period  are  determined  by  a  call  to  RESN.  And,  finally,  the  phase 
portrait  is  plotted  by  making  successive  calls  to  CONTUR  to  obtain  the 
data  points. 

The  actual  plotting  was  done  on  a  CALCOMP  plotter  and  was 
facilitated  by  the  availability  of  an  extensive  library  of  plot  routines 
on  the  CDC  computer. 
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Subroutine  RESN.  RESN  is  the  routine  responsible  for  calculation 


of  many  of  the  parameters  associated  with  the  resonance  term  of 
interest.  It  begins  by  determining  the  nominal  two-body  radius  for  a 
semi-synchronous  orbit.  Then,  using  a  simple  binary  search  algorithm, 
it  searches  out  the  location  of  the  unstable  and  stable  equilibrium 
points  based  on  the  equations  developed  in  Chapter  II.  Next,  again 
using  a  binary  search  algorithm,  it  locates  the  edges  of  the  stable 
resonance  band  and  determines  the  resonance  width.  It  does  this  by 
taking  the  value  of  the  Hamiltonian  at  the  unstable  point  and,  since  the 
energy  contours  passing  through  the  unstable  point  define  the  edges  of 
the  stable  region,  it  searches  along  the  radius  through  the  stable  point 
for  the  same  values.  Finally,  it  determines  the  libration  period  for 
small  oscillations  about  the  stable  point. 

Subroutine  CONTUR.  The  CONTUR  subroutine  is  the  real  heart  of  the 
phase  portrait  plotting.  The  routine  was  provided  by  Dr.  William  Wiesel 
of  the  Astronautics  Department.  Since  the  variables  to  be  plotted  on 
the  phase  portrait,  the  conjugate  momentum  S  and  coordinate  s,  are 
action-angle  variables,  the  routine  did  have  to  be  modified  slightly  to 
permit  contouring  of  a  function  in  polar  coordinates  instead  of 
Cartesian. 

-The  routine  itself  uses  a  two  variable  Newton  iteration  scheme  to 
locate  successive  values  of  S  and  s  for  a  constant  energy  contour  of  the 
Hamiltonian  given  initial  values  of  S  and  s  on  the  contour.  The 
derivatives  or  slopes  of  the  Hamiltonian  are  used  to  determine  the 
direction  the  routine  needs  to  iterate  in  order  to  follow  the  contour. 

Subroutine  FDF.  The  routine  FDF  implements  the  equations  developed 
in  Chapter  II  for  each  resonance  term  of  interest.  Given  values  of  the 
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variables  S  and  s  as  well  as  the  constants,  it  evaluates  the  value  of 


the  coefficients  A  and  B,  the  Hamiltonian,  and  its  derivatives  for  use 
by  the  calling  routine. 

Procedure 

Investigation  of  the  three  resonance  terms  actually  began  with  the 
software  debugging  process  which  continued  throughout  the  study  due  to 
quirks  peculiar  to  each  term.  After  numerous  computer  runs,  though,  a 
fair  amount  of  confidence  in  the  approach  was  developed. 

One  of  the  first  decisions  to  be  made  was  the  choice  of  the 
geopotential  model.  Since  differences  among  models  are  generally  minor, 
the  1973  Smithsonian  Standard  Earth  (Ref  12)  was  chosen  because  of  its 
ready  availability.  Table  I  gives  a  listing  of  the  requisite  constants 
based  on  this  model. 


Table  I 


SAO 

III  Geopotential 

Constants 

Constant 

Symbol 

Value 

Earth  Radius 

R 

e 

6378140  meters 

Canonical  Time  Unit 

TO 

806.8108  sec 

Earth  Rotation  Rate 

«« 

0.0588335721  rad/TU 

First-Order  Zonal  Term 

1082.637  x  10"* 

Second  Order  Zonal  Terra 

-  1.617999  x  10'* 

Principal  Tesseral  Term 

2 

2.7438636  x  10~‘ 

Tesseral  Longitude 

*12 

-  14.923723' 

A  related  question  was  the  values  to  be  assigned  to  the  constants 
of  motion  Q  and  R.  Since  the  second  transformation  relates  Q,  R  and  S 
to  the  orbital  inclination,  a  value  of  zero  was  chosen  for  R  so  as  to 
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result  in  consideration  of  inclinations  very  near  but  not  at  90°  for 
small  values  of  S.  The  exact  value  of  the  inclination  at  an  equilibrium 
point  then  depends  on  Q  and  the  resulting  value  of  S  at  the  point.  The 
choice  of  Q  was  equally  as  arbitrary.  Since  Q  and  S  are  also  related  to 
the  eccentricity,  values  of  Q  for  each  resonance  were  chosen  that  would 
result  in  a  value  of  S  at  the  stable  point  that  corresponded  to  an 
eccentricity  of  about  0.01.  This  was  considered  to  be  both  a  realistic 
and  realizable  eccentricity  for  satellites  in  this  type  of  orbit.  This 
proved  to  be  an  iterative  process  because  of  the  impact  of  the  resonance 
structure  itself  and  the  resultant  plotting  considerations. 

It  was  also  decided  at  this  time  to  set  the  J22  longitude  term, 

X2 2 ,  to  zero.  This  was  done  primarily  as  a  plotting  convenience.  Since 
this  term  only  serves  to  rotate  the  orientation  of  the  phase  portrait  in 
terms  of  the  angle  s,  its  inclusion  adds  little  insight  to  the  problem. 
It  could  just  as  well  have  been  incorporated  in  the  critical  argument  of 
the  second  transformation. 

The  next  step  was  to  run  a  simplified  version  of  RESN  by  itself  in 
order  to  determine  the  approximate  location  of  the  resonance  structures 
prior  to  actually  attempting  to  plot  them. 

With  these  preliminary  steps  taken  care  of,  full  scale 
investigation  of  the  resonance  terms  proceeded.  The  actual  process 
consisted  of  making  "several"  full  scale  computer  runs  on  each  resonance 
term  in  order  to  select  a  Q  and,  hence,  the  resulting  values  of  S  for 
each  term.  The  outputs  of  these  runs  were  some  of  the  resonance 
parameters  associated  with  each  term  as  well  as  polar  contour  plots  of 
the  conjugate  action-angle  variables  S  and  s  for  various  values  of  the 
Hamiltonian. 
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This  process  required  a  certain  amount  of  manual  interaction 
because  of  the  individual  resonance  quirks  mentioned  previously.  For 
example:  the  first  resonance  required  plotting  -S  because  of  the  second 

transformation.  The  second  resonance  was  fairly  routine  but  the  third 
required  major  massaging  because  of  its  extremely  fine  structure. 

The  specific  results  stemming  from  this  process  are  detailed  in  the 


next  chapter. 


IV  Results 


The  specific  results  of  this  investigation  are  presented  in  the 
following  pages.  They  are  grouped  by  resonance  term  as  defined  by  the 
equations  of  Chapter  II.  Tables  II,  III,  and  IV  list  the  parameters 
associated  with  each  resonance  term  and  Figures  1,  2,  and  3  are  the 
corresponding  phase  portraits. 

Note  that  the  values  of  Q  and  S  are  expressed  in  canonical  units 
where  Q+S  =  / ya,  a  in  terms  of  earth  radii  (ER).  Equilibrium  point 
radius  is  also  given  in  ER.  This  is  equivalent  to  the  semi-major  axis 
at  the  equilibrium  point.  The  distance  of  the  equilibrium  points  from 
nominal  is  referenced  to  a  nominal  two-body  semi-synchronous  radius  of 
4.1645031  ER. 

The  phase  portraits  are  standard  polar  contour  plots  of  the  action 

(momentum)  variable  S  (or  -S)  against  the  angle  variable  s  for  equal 

* 

values  of  the  system  Hamiltonian,  F  .  Examination  of  the  plots  indicate 
the  existence  of  one  stable  and  one  unstable  equilibrium  point  (marked 
by  Xs)  associated  with  each  term.  The  stable  point  is  distinguished  by 
the  closed  contours  around  it. 

An  interpretation  and  complete  discussion  of  these  results  and  the 
resonance  effects  is  presented  in  the  next  chapter. 
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Table  II 


First  Resonance  Term 
(Q  =  2.04079386) 

Unstable  Equilibrium  Point 
S  = 
s  = 

Radius  = 

Nominal  Displacement  = 

Stable  Equilibrium  Point 
S  = 

8  * 

Radius  = 

Nominal  Displacement  = 

Orbit  Inclination  (I)  = 

Orbit  Eccentricity  (e)  = 

Resonance  Width  » 

Libration  Period  = 


0.00007779  ERVTU 
180.0° 
4.16452206  ER 
120.8145  meters 

0.00008679  ER*/TU 
0.0° 

4.16448533  ER 
-  113.4052  meters 
90.005° 
0.009 

2791.9897  meters 
12,188.8  days 
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1ST  RESONANCE  T  F  R 


M 


( Q=2. 04079366 


Figure  1.  First  Resonance  Term 


Table  III 


Second  Resonance  Tern 

(Q  =  2.04048689) 

Unstable  Equilibrium  Point 
S  = 
s  = 

Radius  = 

Nominal  Displacement  = 

Stable  Equilibrium  Point 
S  * 
s  * 

Radius  = 

Nominal  Displacement  * 

Orbit  Inclination  (I)  = 

Orbit  Eccentricity  (e)  » 

Resonance  Width  ° 

Libration  Period  = 


0.00020756  ER’/TU 
0.0° 

4.16443385  ER 
-  44’l.7858  meters 

0.00024027  ER*/TU 
180.0° 
4.16456734  ER 
409.6754  meters 
89.987° 
0.015 

8740.8379  meters 
3,843.1  days 
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2ND  RESONANCE  TERM 

( Q-2. 04043689  1 

Figure  2.  Second  Resonance  Terra 
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Table  IV 


Third  Resonance  Term 
(Q  =  2.04068145) 

Unstable  Equilibrium  Point 
S  = 

s  = 

Radius  = 

Nominal  Displacement  = 

Stable  Equilibrium  Point 
S  = 

8  = 

Radius  = 

Nominal  Displacement  * 

Orbit  Inclination  (I)  = 

Orbit  Eccentricity  (e)  = 

Resonance  Width  = 

Libration  Period  = 


0.0002996  ERl/TU 
0.0° 

4.16450307  ER 

-  0.2783  meters 

0.00002996  ER*/TU 
180.0° 
4.16450307  ER 

-  0.2708  meters 

89.998° 
0.009 
5.546  meters 
6,385,139.3  days 
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V  Discussion  and  Recommendations 


Discussion 

One  of  the  goals  of  this  study  was  to  determine  if  consideration  of 
all  resonance  terms  associated  with  a  particular  tesseral  harmonic  and 
commensurabi 1 ity  ratio  would  result  in  multiple  usable  equilibrium 
points.  From  a  first  glance,  the  results  do  not  look  particularly 
promising.  Investigation  of  the  structure  of  all  three  resonance  terms 
indicate  the  existence  of  only  a  single  stable  point  associated  with 
each  for  non-circular  orbits.  In  addition,  two  other  aspects  of  the 
stable  points  raise  questions  about  their  suitability. 

First,  note  that  all  three  stable  points  are  displaced  slightly 
from  nominal  semi-synchronous  altitude.  The  broader  and  stronger  the 
resonance,  the  greater  the  displacement.  This  results  in  objects  placed 
at  each  stable  point  having  mean  motions  different  from  each  other  as 
well  as  from  nominal  semi-synchronous  causing  potential  phase  drifts 
between  the  different  terms.  The  actual  physical  implications  of  this 
will  be  made  clear  later. 

Secondly,  the  stable  region  due  to  the  third  term  is  so  shallow  and 
narrow  that  it  may  prove  impossible  to  use.  With  the  resonance  being 
only  5.5  meters  wide,  it  would  be  difficult  to  position  a  satellite  in 
that  small  a  region  at  that  altitude  to  begin  with.  Even  if  one  were 
placed  there,  the  resonance  is  so  weak,  being  of  O(e’),  that  it  would 
require  only  a  minor  external  perturbation  to  push  it  out,  thereby 
negating  its  value. 

Although  the  third  term  does  not  appear  a  likely  candidate  for  use, 
a  closer  look  at  the  first  and  second  terms  may  be  merited.  This  is 
because  that,  while  the  position  of  a  stable  point  in  terms  of  the 
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momentum  variable  S  is  easily  interpreted  as  determining  the  semi-major 
axis  of  the  resonant  orbit,  the  physical  meaning  of  the  angle  variable  s 
is  not  as  clear-cut.  Because  of  the  second  transformation,  s  actually 
incorporates  several  different  non-parallel  orbital  angles,  each  of 
which  may  possess  secular  variations.  This  makes  it  difficult  to 
visualize  just  what  the  resonant  condition  is  for  a  polar  orbit  without 
closer  study. 

First  Resonance  Term.  Returning  to  Eq  (13)  to  express  s  in  terms 
of  the  Delaunay  elements  gives 

tj  +  2h ,  -  2n,t  +  2g,  -  2X j 2  =  0 

for  the  stable  point  of  the  first  resonance.  Initially  setting  g,  =  0 
and  solving  for  the  longitude  of  the  node  by  setting  li  =  0  gives 

hi  -  n0t  =  Xj t 

where  h  -  n,t  is  the  geocentric  longitude  of  the  ascending  node.  In  the 
same  manner,  solving  for  the  longitude  of  the  descending  node  by  setting 
tj  =  *  gives 


In  other  words,  resonant  stability  for  this  case  is  determined  not  only 
by  the  semi-major  axis  of  the  orbit  but  by  also  locating  the  ascending 
node  of  the  orbit  over  the  major  axis  of  the  earth's  equatorial  section 
and  the  descending  node  over  the  minor  axis.  This  is  certainly 
intuitively  satisfying  from  a  physical  point  of  view.  A  similar 
analysis  performed  on  the  unstable  point  where 

l,  +  2h,  -  2n,t  +  2g,  -  2X,  2  ■  ir 

^  gives 

i 

4* 
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h,  -  n,t  =  X,,  ♦  j 

for  the  longitude  of  the  ascending  node  and 

hj  -  n0 1  =  X2  2 

for  the  descending  node.  This  indicates  that  unstable  resonance  results 
from  the  first  term  when  the  ascending  node  of  the  orbit  is  over  the 
minor  axis  and  the  descending  node  is  over  the  major  axis.  Unlike  the 
geosynchronous  equatorial  case  where  resonance  is  a  function  of  semi¬ 
major  axis  and  satellite  longitude,  in  the  polar  case  it  is  a  function 
of  semi-major  axis  and  node  location. 

The  preceding  analysis  was  based  on  an  argument  of  perigee  of  zero. 
Now  consider  what  happens  if  some  other  argument  of  perigee  is  selected. 
If  gi  =  ir/4  then 

l,  +  2h ,  -  2n0 t  +  |  -  2X22  =  0 
for  the  stable  point.  This  yields 

h,  -n,t=  X  j  2  -  ^ 

for  the  longitude  of  the  ascending  node  and 

,  *  -  3* 

hl-n,t-A2J-^ 

for  the  longitude  of  the  descending  node.  Again,  this  is  understandable 
from  a  physical  point  of  view.  Since  a  rotation  of  the  perigee  causes 
the  time  from  ascending  node  to  descending  node  to  differ  from  the  time 
from  descending  to  ascending,  a  change  in  the  node  is  necessary  to 
balance  the  effects  of  the  earth's  equatorial  section  and  maintain  the 
resonance  condition.  In  general,  a  rotation  of  the  perigee  location 
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results  in  an  equivalent  regression  of  the  line  of  nodes  such  that,  for 

gi  =  */2 

it 

h,  -n0t-  XJ2  2 

for  the  longitude  of  the  ascending  node  and 

h i  —  n0 1  =  X  j  2  ~  * 

for  the  descending  node  so  that  resonant  stability  is  now  determined  by 
an  ascending  node  over  the  longitude  of  the  minor  axis  and  descending 
node  over  the  major  axis.  A  further  rotation  of  perigee  to  ir  results  in 
ascending  and  descending  nodes  being  on  the  opposite  ends  of  the  major 
and  minor  axes  then  for  the  g,  =  0  case.  A  similar  analysis  on  the 
unstable  equilibrium  point  shows  the  same  effect. 

This  opens  up  an  array  of  opportunities  for  deploying  satellite 
systems.  Simply  by  careful  selection  of  the  longitude  of  ascending  node 
and  the  argument  of  perigee  for  each  orbit,  a  constellation  of  almost 
any  size,  depending  on  mission  requirements,  can  be  established  that 
utilizes  the  first  term  resonance  effects.  Orbits  of  common  nodal 
longitude  can  be  used  to  provide  the  same  ground  tracks  while  different 
arguments  of  perigee  will  determine  the  particular  node  and,  hence,  its 
ground  track.  In  addition,  since  they  all  depend  on  the  first  term, 
they 'all  have  the  same  semi-major  axis  and,  as  a  result,  the  same  mean 
motion,  resulting  in  a  constant  phasing  among  the  satellites  of  the 
constellation.  The  last  issue  raises  one  final  consideration.  The 
previous  analysis  has  shown  that  resonance  for  a  semi-synchronous  near 
polar  orbit  is  fundamentally  a  locking  of  the  nodes  of  the  orbit  in 
phase  with  the  axes  of  the  earth's  equatorial  section.  It  was  pointed 
out  earlier,  however,  that  the  altitudes  of  the  stable  points  are 


displaced  slightly  from  nominal  semi-synchronous  altitude,  slightly 
below  in  the  case  of  the  first  resonance.  An  object  placed  at  the 
stable  semi-major  axis  of  the  first  resonance  will  have  a  mean  motion  of 
2.00001281  revolut ions/day .  This  will  result  in  a  difference  between 
the  longitude  of  the  nodes  and  the  major  and  minor  axes  of  0.0023°/day. 
As  a  result,  one  effect  of  the  resonance  is  to  cause  a  secular  variation 
in  the  right  ascension  of  -0.0023°/day  in  order  to  maintain  the  phasing 
of  the  nodes  with  the  axes  of  the  earth's  equator.  In  truth,  because  of 
the  interaction  of  h  and  g  for  the  first  term,  this  variation  is  imposed 
on  both  as  a  drift  in  both  h  and  g,  the  exact  components  of  which  are 
more  difficult  to  separate.  This  effect,  however,  provides  an 
explanation  of  why  the  stronger  the  resonance  and,  hence,  the  greater 
the  secular  variation,  the  greater  the  displacement  from  nominal 
altitude. 

Second  Resonance  Term.  Taking  the  argument  of  the  second  term  from 
Eq  (13)  gives 

l*  +  2hj  -  2n0t  -  2\tt  =  ir 
for  the  stable  point.  Setting  Z*  =  0  then  yields 


as  the  longitude  of  the  ascending  node.  In  the  same  manner,  letting 
Zj  *  w  results  in 

hj  —  n,  t  *  Xj2 

for  the  longitude  of  the  descending  node.  Therefore,  for  the  second 
resonance  term,  resonant  stability  is  established  by  placing  ascending 
node  over  the  minor  axis  and  descending  node  over  the  major.  Performing 
the  same  analysis  on  the  unstable  point  where 
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2n0t  -  2X j  2 


0 


gives 


*1  +  2h2  - 


for  ascending  node  and 


h2  n8 1  a 


hj  -  n„t  = 


ir 

2 


for  descending  node  or  ascending  at  major  axis  and  descending  at  minor. 
For  the  second  term,  however,  since  the  argument  of  perigee  does  not 
enter  into  the  argument,  rotation  of  the  perigee  does  not  alter  the 
resonant  condition.  As  a  result,  stable  resonance  is  always 
characterized  by  the  fixed  relationship  of  the  nodes  to  the  equatorial 
axes.  This  makes  the  second  term  somewhat  less  interesting  than  the 
first. 

Now,  the  stable  altitude  of  the  second  term  also  imposes  a  secular 
variation  on  the  line  of  nodes.  Its  semi-major  axis  is  equivalent  to  a 
mean  motion  of  1.99995373  revolut ions/day .  As  a  result,  the  right 
ascension  must  precess  0.0083°/day  in  order  to  maintain  the  resonant 
condition.  In  this  case,  the  effect  is  totally  in  right  ascension  since 
argument  of  perigee  does  not  appear  in  the  argument. 

Third  Resonance  Term.  For  completeness,  the  third  term  will  be 
brief-ly  addressed.  Again,  taking  the  argument  from  Eq  (13)  for  the 
stable  point  gives 

£i  +  2h,  -  2n0t  -  2g,  -  2A2j  =  ir 
Setting  g,  ■  0  results  in 


for  the  longitude  of  ascending  node  and 
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for  descending  node.  As  a  result,  the  basic  resonant  condition  is 
identical  to  that  of  the  second  term,  however,  since  the  argument  of 
perigee  appears  in  the  argument,  a  rotation  of  perigee  results  in  a 
rotation  of  the  line  of  nodes  similar  to  that  of  the  first  term  but  in 
the  opposite  direction.  The  stable  altitude  of  the  third  point  gives  a 
mean  motion  of  2.00000003  revolutions/day  which  results  in  a  resonant 
variation  in  right  ascension  of  0.000005 6° /day  in  order  to  stay  in 
phase. 

Because  of  its  closer  proximity  to  a  nominal  semi-synchronous 
altitude,  the  third  term  might  be  the  most  interesting  were  it  not  for 
its  relative  weakness.  In  fact,  because  of  the  widths  of  the  first  and 
second  resonances  which  would  encompass  the  third  for  a  common  resonant 
longi-tude,  an  object  in  the  third  would  more  likely  be  dominated  by  the 
stronger  term  and  tend  to  librate  within  that  resonance.  The  “<ame 
result  could  occur  for  the  first  term,  since  the  8700  meter  width  of  the 
second  term  totally  contains  the  first  at  a  common  stable  longitude. 

Libration  Periods .  This  last  thought  raises  the  point  of  libration 
again.  From  the  numerical  results  presented  in  the  last  chapter,  it  can 
be  seen  that  all  three  resonance  terms  have  relatively  long  libration 
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periods  making  the  orbital  changes  due  to  libration  small  over  short 
periods  of  time.  Because  of  the  interpretation  just  given  for  the 
resonance  variables,  the  libration  is  actually  seen  as  an  oscillation  in 
the  semi-major  axis  and  longitude  of  ascending  node  of  the  orbit  for  a 
satellite  not  placed  exactly  at  a  stable  point.  The  further  the 
displacement  in  either  direction  from  the  stable  point,  the  longer  the 
libration  period  until  it  falls  outside  the  stable  region  where  it 
exhibits  a  secular  motion  in  longitude  known  as  circulation. 

Although  not  a  lot  of  data  is  available  on  semi-synchronous 
libration  periods,  a  couple  of  studies  have  obtained  analytical  values 
for  other  inclinations  that  provide  some  comparison.  Allan  (Ref  2:63) 
obtained  a  value  of  approximately  2265  days  for  a  semi-synchronous  orbit 
with  an  inclination  of  34.4°.  Sochilina  (Ref  20:345-348)  gives 
approximately  5000  days  as  the  librationai  period  at  the  critical 
inclination.  These  appear  to  compare  favorably  with  the  value  of  3851 
days  for  the  dominant  second  term  obtained  numerically  in  this  study. 

Recommendations 

From  the  preceding,  it  is  obvious  that  there  does  exist  a  definite 
opportunity  for  exploiting  resonance  effects  on  a  semi-synchronous  near 
polar  orbit.  As  was  shown  in  the  discussion,  the  effect  of  the  argument 
of  perigee  in  the  first  resonance  term  provides  significant  flexibility 
in  the  deployment  of  satellites  utilizing  resonance  effects.  This  is 
something  that  has  been  overlooked  in  the  past  because  of  the  tendency 
to  consider  only  circular  orbits,  thus  eliminating  any  resonance  effects 
at  all,  or  to  treat  only  the  dominant  term  which  in  this  case  does  not 

f  contain  the  argument  of  perigee  effect.  As  a  result,  this  study  has 

m  - 

shown  that  consideration  of  all  resonance  terms  does  provide  additional 


42 


stable  positions  which  because  of  their  physical  interpretation  do  not 
interfere  with  each  other  as  much  as  might  be  expected.  These  results 
certainly  merit  consideration  in  the  deployment  of  future  satellite 
systems.  In  this  regard,  there  are  two  additional  areas  worth  looking 
at . 

One  would  be  to  study  the  added  effects  of  additional  *-esseral 
harmonics.  Even  though  they  are  substantially  smaller  that  the  J22 
term,  they  may  affect  the  results  enough  to  require  consideration. 

Another  area  of  investigation  that  this  study  had  hoped  to  address 
was  to  determine  the  station  keeping  requirements  in  the  vicinity  of  the 
stable  points.  This  could  prove  to  be  the  real  test  of  the  feasibility 
of  using  these  regions  for  the  deployment  of  future  satellite  systems. 
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Appendix 

Program  and  Subroutine  Listings 

The  following  is  a  sample  listing  of  the  program  and  three  main 
subroutines  (set  up  for  the  first  resonance)  used  in  the  process  of  this 
investigation.  The  routines  called  for  plotting  are  from  a  general 
purpose  plotter  subroutine  library  for  the  CDC. 


PROGRAM  RESPLT 

Q  ***************** 

C  DRIVER  PROGRAM  FOR  INVESTIGATION  AND  PLOTTING 
C  OF  RESONANCE  EFFECTS  DUE  TO  EARTH'S  GEOPOTENTIAL 

C  ***************** 

DIMENSION  BIGS(9000),  SMLS(9000) 

REAL  J2, J2SQ, J4, J40VR, J22.N0 
COMMON  /SAOIII/  J2, J2SQ, J4, J40VR, J22.N0, PHI22 
COMMON  /PARAM/  Q . R, RADIAN , TP I 
C  SET  GEOPOTENTIAL  AND  HAMILTONIAN  CONSTANTS 
J2  =  1082.637E-6 
J2SQ  *  J2**2 
J4  *  -1.617999E-6 
J40VR  =  J4/J2SQ 
J22  =  2.7438636E-6 
NO  =  5.86729371E-2 
PHI22  =  0.0 
Q  =  2.04265452 
R  =  0.0 

PI  *  4.0*ATAN(1.0) 

TP I  =  2.0  *  PI 
RADIAN  =  360. /TPI 
C  SCALE  AND  INITIALIZE  PLOT 
SCALE  *  .00005 
CALL  PLOTS (0,0 , 1 ) 

CALL  PLOT(0.0,-.5,-3) 

•  CALL  PLOT(0. 0,0. 0,-3) 

CALL  PL0T( 1 .0, 7.0, 3) 

CALL  PLOT(4 .0 , 7.0,2) 

CALL  AXIS(4. 0,7.0, ' (-S) ' ,-4, 3. ,0. ,0.0, SCALE) 

CALL  PLOT(4 .0 ,4 .0 , 3) 

CALL  PLOT(4. 0,10. 0,2) 

CALL  PLOT( 1 .2 ,2 . 5 , 3) 

CALL  SYMBOL( 1,23,2.5,0.315, '1ST  RESONANCE  TERM', 0., 18) 
CALL  SYMBOL(3. 05, 2. 0,0. 14, ' (Q=2. 04265452)'  ,0., 14) 

CALL  PLOT(0.0,0.0, 3) 

SI  *  -.000025 
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C  CALCULATE  AND  PLOT  EQUILIBRIUM  POINTS 

CALL  RESN(SI,STBLE,UNSTB) 

XS  =  4.0  +  ABS(STBLE) /SCALE 
XU  =  4.0  -  ABS(UNSTB) /SCALE 
CALL  SYMBOL(XS ,7.0,0.07,4,0. ,-l ) 

CALL  SYMBOL ( XU ,7 .0,0.07,4,0. ,-l  ) 

0  ***************** 

C  LOOP  TO  PLOT  CONTOURS 

C  ***************** 

DO  4  I  =  1,10 

C  INITIALIZE  FOR  INDIVIDUAL  CONTOUR 

BIGS ( 1 )  =  -SI 
BSNEW  =  SI 
SMLS(l)  =  0.0 
SMLSNW  =0.0 
IFLG  =  0 
N  =  1 
ICAS  =  0 
DEL  =  TPI/3600. 

C  DETERMINE  POINTS  ON  EACH  CONTOUR 

1  CALL  CONTUR( BSNEW .SMLSNW, DEL, BSNEW, SMLSNW, IOK, ICAS, FS) 

IF  (IOK.EQ.O)  GO  TO  2 

PRINT  *, 'BGSO=' ,BIGS(N), '  BIGS= ', BSNEW, '  SMLS=’ .SMLSNW, '  FSTR=',FS 
GO  TO  3 

2  N  =  N  +  1 
BIGS(N)  =  -BSNEW 
SMLS(N)  =  SMLSNW  *  RADIAN 
IF  (N.EQ.8998)  GO  TO  3 

C  END  OF  CONTOUR? 

IF  ((SMLSNW. GE. 0.0). AND. (IFLG. EQ.l))  GO  TO  3 
IF  (SMLSNW. LT. 0.0)  IFLG  =  1 
IF  (SMLSNW. LT.TPI)  GO  TO  1 

3  BIGS(N+1 )  =  0.0 
BIGS(N+2)  =  SCALE 

Q  ***************** 

C  PLOT  EACH  CONTOUR 

Q  ***************** 

CALL  POLAR( BIGS ,SMLS ,-N,4 , 0,7. 0,3. 0,0, 3) 

C  SET  BIG  S  FOR  START  OF  NEXT  CONTOUR 

SI  =  SI  -  .00001 
.  IF  (I.EQ.7)  SI  =  -.000033 
IF  (I.EQ.8)  SI  *  -.000141 
IF  (I.EQ.9)  SI  «  -.00015 

4  CONTINUE 
K  «  0 

C  TERMINATE  PLOT 

CALL  PLOTE(K) 

STOP 

END 
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SUBROUTINE  RESN(St ,STBLE , UNSTB) 

***************** 

ROUTINE  TO  CALCULATE  RESONANCE  PARAMETERS 
SI  IS  INITIAL  BIG  S  TO  START  SEARCH 
***************** 

REAL  J2, J2SQ, J4, J40VR, J22.N0 
COMMON  /SAOIII/  J2, J2SQ, J4, J40VR, J22,N0,PHI22 
COMMON  /PARAM/  Q , R , RADIAN , TP  I 
C  INITIALIZE  CONSTANTS  AND  VARIABLES 

ER  =  6378140.0 
BGS  =  SI 
DS  =  -7  *  SI 
SMS  =  TP 1/2.0 

C  CALCULATE  NOMINAL  SEMI-SYNCHRONOUS  RADIUS 

SR  =  (2*N0)**(-2. 0/3.0) 

PRINT  *, '  ' 

PRINT  *, 'NOMINAL  RADIUS  (DU)' 

PRINT  *,  SR 

c  ***************** 

C  LOCATE  UNSTABLE  POINT 

C  ***************** 

1  DS  =  DS/2.0 

CALL  FDF( BGS , SMS , FSTR, DFDBS , DFDSS ) 

BGS  =  BGS  -  DS  *  SIGN( 1.0, DFDBS) 

IF  (DS.GT.1E-16)  GO  TO  I 

FUNSTB  =  FSTR 

UNSTB  »  BGS 

RAD  =  (Q  +  BGS)**2 

RADMET  =  (RAD  -  SR)*ER 

PRINT  *,  '  * 

PRINT  *, 'UNSTABLE  EQUILIBRIUM  POINT' 

PRINT  *,  '  ' 

PRINT  *, 'S=' .UNSTB, '  RADIUS  (DU)-', RAD,'  DIST  FM  NOM  (M)-’ .RADMET, 
2  '  FSTR-' .FUNSTB 
C  REINITIALIZE 

BGS  -  SI 
DS  -  -7  *  SI 
SMS  -  0.0 

Q  ***************** 

C  LOCATE  STABLE  POINT 

0  ***************** 

2  DS  -  DS/2.0 

CALL  FDF( BGS , SMS , FSTR , DFDBS , DFDSS ) 

BGS  -  3S  -  DS  *  SIGNQ.O, DFDBS) 

IF  (DS.GT. IE-16)  GO  TO  2 
S^BLE  -  BGS 
RAD  «  (Q  +  BGS )**2 
RADMET  -  (RAD-SR)*ER 
PRINT  *, '  ' 

PRINT  *, 'STABLE  EQUILIBRIUM  POINT' 

PRINT  *, '  ' 

PRINT  *, ' S-' ,STBLE , '  RADIUS  (DU)-', RAD,'  DIST  FM  NOM  (M)-' .RADMET, 
2  '  FSTR-' , FSTR, '  DFDBS- ' , DFDBS , '  DFDSS- ', DFDSS 
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C  ***************** 

CALCULATE  RESONANCE  WIDTH 
***************** 

DS  =  5  *  SI 

LOCATE  INSIDE  EDGE  OF  STABLE  REGION 

3  DS  =  DS/2.0 

CALL  FDF( BGS , SMS , FSTR, DFDBS , DFDSS ) 

FDIF  =  FUNSTB  -  FSTR 
BGS  =  BGS  -  DS  *  SIGN( 1.0, FDIF) 

IF  (ABS(DS) .GT. IE-16)  GO  TO  3 
EDGEI  =  BGS 
PRINT  *,  '  ' 

PRINT  *, 'INSIDE  EDGE*', EDGEI 
ED1R  =  (Q+BGS)**2 
BGS  =  STBLE 
DS  =  5  *  SI 

LOCATE  OUTSIDE  EDGE  OF  STABLE  REGION 

4  DS  =  DS/2.0 

CALL  FDF( BGS , SMS , FSTR, DFDBS , DFDSS) 

FDIF  *  FUNSTB  -  FSTR 
BGS  =  BGS  +  DS  •  SIGN( 1.0, FDIF) 

IF  (ABS(DS).GT. 1C-16)  GO  TO  4 
EDGE 2  =  BGS 
PRINT  *, '  ' 

PRINT  *, 'OUTSIDE  EDGE* ' ,EDGE2 
ED2R  *  (Q  +  BGS )**2 
RW  =  (ED1R  -  ED2R)  *  ER 
PRINT  '  ' 

PRINT  *, 'RESONANCE  WIDTH  (M)=',RW 
***************** 

CALCULATE  LIBRATION  PERIOD 
***************** 

CALCULATE  SMALL  DISPLACEMENT  NEAR  STABLE  POINT 
DELSS  =  TP 1/ 1000. 

DELBS  =  (EDGE2  -  EDGEI )/1000. 

C  NUMERICALLY  DETERMINE  SECOND  PARTIALS 

CALL  FDF( STBLE, SMS, FSTR, S DFDBS .SDFDSS) 

BGS  =  STBLE  +  DELBS 

CALL  FDF( BGS , SMS , FSTR, DFDBS , DFDSS ) 

DFDBS2  »  (DFDBS  -  SDFDBS)/DELBS 
•  SMS  =  SMS  +  DELSS 
CALL  FDF(STBLE, SMS, FSTR, DFDBS, DFDSS) 

DFDSS2  *  (DFDSS  -  SDFDSS)/DELSS 
C  DETERMINE  FREQUENCY  SQUARED 

FREQSQ  *  ABS( DFDBS2  *  DFDSS2) 

C  CALCULATE  PERIOD 

PER  -  N0/SQRT( FREQSQ) 

PRINT  ’ 

PRINT  *, 'DELBS*' , DELBS, '  DELSS*' , DELSS, ’  DFDBS2- ’ .DFDBS2, 
2  *  DFDSS2-' .DFDSS2 
PRINT  *, '  ' 

PRINT  *, 'LIBRATION  PERIOD  (DAYS)*', PER 

RETURN 

END 
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SUBROUTINE  CONTUR(X, Y, DEL ,XNEW, YNEW , IOK, ICAS , FF ) 

C  ***************** 

C  ROUTINE  TO  DETERMINE  NEW  POINT  ON  CONTOUR  GIVEN  OLD 
C  ON  INPUT  X , Y  IS  PT  ON  CONTOUR 

C  FUNCTION  VALUE  ON  CONTOUR  FF  IS  DEFINED  ON  FIRST  CALL 

C  ICAS  IS  ZERO  ON  FIRST  CALL,  USED  TO  FOLLOW  CONTOUR  IN  SAME 

C  SENSE  THROUGHOUT  AS  SLOPE  CASE  SWITCHES 
C  EACH  CALL  UPDATES  TO  NEW  POINT  ON  CONTOUR,  ONE  PER  CALL 

C  ***************** 

TOLERANCES . 

TEN  FIGURES  IN  X  OR  Y  . 

TOLXY  «  IE-10 

.  OR  FOURTEEN  IN  THE  FUNCTION 

TOLF  =  IE-14 
ASSUME  SUCCESS 
IOK  =  0 

FUNCTION  VALUE  AND  PARTIALS 
CALL  FDF(X,Y,F,DFDX,DFDY) 

STORE  FUNCTION  VALUE  ON  FIRST  CALL 
IF(ICAS  .EQ.  0)  FF  =  F 

BRANCH  ON  CASE 

IF(ABS( DFDY/X )  .LT.  ABS(DFDX))  GO  TO  500 
***************** 

SHALLOW  SLOPE,  FIX  X  DETERMINE  Y 
***************** 

FIRST  CALL  ? 

IF(  ICAS  .EQ.  0)  ICAS  »  -1 
HAS  CASE  SWITCHED? 

IF( ICAS  .EQ.  -1)  GO  TO  50 

CASE  SWITCHED  -  FIX  DEL  TO  CONTINUE  TO  FOLLOW  IN  SAME  DIRECTION 
DEL  =  DEL  *  SIGNU.O.DFDX)  *  SIGN(1 .0,DFDY)*(-1 .0) 

ICAS  =  -1 
50  CONTINUE 

DELX  =  DEL  *  ABS(X) 

XNEW  =  X  +  DELX 
YNEW  =  Y 

DY  =  -DFDX*DELX/DFDY 
NEWTON  LOOP 
•  DO  100  I  =  1,30 
YNEW  -  YNEW  +  DY 
RESIDUAL  AND  SLOPE 
CALL  FDF ( XNE W , YNEW , FP , DF  P  DX , DF  P  DY ) 

RES  -  FP  -  FF 
DY  -  -RES/DFPDY 
CONVERGED? 

ABSOLUTE  AND  RELATIVE  FUNCTION  ERROR 
ERRF1  -  ABS(RES) 

ERRF2  -  1.0 

IF(FF  .NE.  0.0)  ERRF2  -  ABS(ERRF1/FF) 

ERRF  -  AMINI(ERRF1 ,ERRF2) 

IF(ERRF  .LT.  TOLF)  GO  TO  1000 
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C  ABSOLUTE  AND  RELATIVE  POSITION  ERROR 

ERRXYI  =  ABS(DY) 

ERRXY2  =  1.0 

IF( YNEW  .NE.  0.0)  ERRXY2  =  ABS( DY/YNEW) 

ERRXY  =  AMIN1 (ERRXY1 ,ERRXY2 ) 

IF(ERRXY  .LT.  TOLXY )  GO  TO  1000 
100  CONTINUE 
C  FAILURE 

I OK  =  -1 
RETURN 

***************** 

STEEP  SLOPE  CASE,  FIX  Y,  DETERMINE  X 
***************** 

FIRST  CALL  ? 

500  IF(ICAS  .EQ.  0)  ICAS  =  1 
SWITCHED? 

IF( ICAS  .EQ.  1)  GO  TO  501 
IT  SWITCHED 

DEL  =  DEL  *  SIGN(l.O.DFDX)  *  SIGN( 1 .0 , DFDY)*(-1 .0 ) 
ICAS  =  1 

501  CONTINUE 
XNEW  =  X 
YNEW  =  Y  +  DEL 
DX  =  -DFDY*DEL/DFDX 
DO  600  I  =  1,30 
XNEW  -  XNEW  +  DX 

CALL  FDF( XNEW, YNEW, FP.DFPDX.DFPDY) 

RES  =  FP-FF 
DX  =  -RES/DFPDX 
CONVERGED  ? 

ABSOLUTE  AND  RELATIVE  FUNCTION  ERROR 
ERRF1  =  ABS(RES) 

ERRF2  =1.0 

IF(FF  .NE.  0.0)  ERRF2  =  ABS(ERRF1/FF) 

ERRF  =  AMIN1 (ERRF1 ,ERRF2) 

IF(ERRF  .LT.  TOLF)  GO  TO  1000 
C  ABSOLUTE  AND  RELATIVE  POSITION  ERROR 
ERRXYI  *  ABS(DX) 

ERRXY2  =  1.0 

IF (XNEW  .NE.  0.0)  ERRXY2  =  ABS(DX/XNEW) 

•  ERRXY  =  AMIN 1 (ERRXYI ,ERRXY2) 

IF(ERRXY  .LT.  TOLXY)  GO  TO  1000 
600  CONTINUE 
C  FAILURE 
IOK  =  -1 
RETURN 

C  SUCCESS 
1000  CONTINUE 
RETURN 
END 
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SUBROUTINE  FDF( BGS ,SMS , FSTR, DFDBS , DFDSS) 

***************** 

ROUTINE  TO  EVALUATE  HAMILTONIAN  AND  FIRST  PARTIALS 
BGS  AND  SMS  ARE  INPUT  CONJUGATE  ACTION-ANGLE  VARIABLES 
***************** 

REAL  J2.J2SQ, J4, J40VR, J22.N0 
COMMON  /SAOIII/  J2, J2SQ, J4, J40VR, J22.N0, PHI22 
COMMON  /PARAM/  Q, R, RADIAN, TP  I 
QPS  =  Q+BGS 

COMPUTE  COEFFICIENTS  FOR  EQUATIONS 
A  =  (Q+BGS)/ (Q+2*BGS) 

B  =  (2*BGS+R)/ (Q+2*BGS) 

***************** 

FIRST  RESONANCE  TERM  HAMILTONIAN 
C  ***************** 

FSTR  =  l/(2*QPS**2)+J2*A**3/(2*QPS**6)*(3/2*B**2-l/2) 

2  +J2SQ*A**5/(4*QPS**10)* 

3  (15/32+27/32*J40VR-(27/16+135/16*J40VR)*B**2+(I5/32+315/32*J40VR) 

4  *B**4+3/ 8*A* ( 1-6*B**2+9*B**4 )-A**2* (15/ 32+45/ 32* J40VR 

5  -(15/16+225/ 1 6*J40VR)*B**2- (105/ 32-525/ 32*J40VR)*B**4 ) ) 

6  -3*J22*SQRT(-BGS)/(4*SQRT(2.0)*QPS**6.5)*(1+B)**2* 

7  COS(SMS-2*PHI22)+2*NO*BGS 

Q  ***************** 

C  PARTIAL  OF  F*  WITH  RESPECT  TO  BIG  S 
C  ***************** 

DFDBS  =  -l/QPS**3-3/2*J2*A**3/QPS**7*((3*Q+4*BGS)/(Q+2*BGS)* 

2  (3/ 2VB**2-1/ 2)-A*B*(2*Q-2*R)/ (Q+2*BGS) )-5/ 4* J2SQ*A**5/ 

3  QPS**1 1*( ( 3*Q+4*BGS )/ (Q+2*BGS )* 

4  ( 15/32+27/32*J40VR-(27/ 16+1 35/ 16*J40VR)*B**2+( 15/ 32+315/ 32* J40VR) 

5  *B**4+3/ 8*A* ( 1-6*B**2+9*B**4 )-A**2* (15/ 32+45/ 32* J40VR 

6  -( 15/16+225/16* J40VR)*B**2- ( 105/32-525/ 32*J40VR)*B**4 ) ) 

7  -A/5*( (2*Q-2*R)/ (Q+2*BGS)*(-(27/8+l 35/ 8* J40VR)*B+( 1 5/ 8+315/ 8* 

8  J40VR)*B**3+3/ 8*A*(36*B**3-1 2*B)+A**2*( ( 15/ 8+225/ 8*J40VR)*B 

9  +(105/8-525/8*J4OVR)*B**3))-Q/(Q+2*BGS)*(3/8*(l-6*B**2+9*B**4) 

A  -2*A*( 15/32+45/ 32*J40VR-( 15/ 16+225/ 16* J40VR)*B**2 

B  - ( 105/ 32-525/ 32* J40VR)*B**4 ) ) ) ) 

C  -3*J22/ (8*SQRT( 2*(-BGS ) )*QPS**7 . 5 )*( ( 12*BGS-Q )*( 1+B)**2 
D  -(2*Q-2*R)/(Q+2*BGS)*4*BGS*A*(1+B) )*C0S(SMS-2*PHI22)+2*N0 

(J  ***************** 

C  PARTIAL  OF  F*  WITH  RESPECT  TO  SMALL  S 
C  ***************** 

DFDSS  »  3*J22*SQRT(-BGS)/(4*SQRT(2.0)*QPS**6.5)* 

2  ( 1+B)**2*S IN(SMS-2*PHI22) 

RETURN 

END 
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